Introduction

Les fichiers sont déposées ici

l’ensemble du travail est sur le serveur de l’IFB: shared/projects/dubii2021/lnegroni/EvaluationM4M5-main

Travail sur des données de reséquençage d’un génome bactérien : Bacillus subtilis. Les données sont issues de cet article :

Analyses

Organisation de votre espace de travail

# connexion a partir de son PC à IFB
ssh -X lnegroni@core.cluster.france-bioinformatique.fr
#1- Aller dans le bon directory de l'IFB 
cd /shared/projects/dubii2021/lnegroni
# telecharger le TP
wget https://github.com/DU-Bii/EvaluationM4M5/archive/refs/heads/main.zip
# decompresser le .zip
unzip main.zip
#aller dans le dossier cree
cd ./EvaluationM4M5-main
#creer les dossiers de travail
mkdir ./FASQ ./CLEANING ./QC ./MAPPING

Téléchargement des données brutes

Récupérez les fichiers FASTQ issus du run SRR10390685 grâce à l’outil sra-tools [1]

# charger le module de sra-tools
module load  sra-tools
#noter la version
fasterq-dump -h #"fasterq-dump" version 2.10.3
#telecharger les .fasq dans le dossier FASQ
srun --cpus-per-task=1 fasterq-dump --split-files -p SRR10390685 --outdir FASQ
# compresser les fichiers
srun gzip *fasq

Combien de reads sont présents dans les fichiers R1 et R2 ?

#utilisation du module seqkit 
module load seqkit
srun seqkit stats --threads 1 *.fastq.gz
#resultats:
#file                    format  type   num_seqs        sum_len  min_len  avg_len  max_len
#SRR10390685_1.fastq.gz  FASTQ   DNA   7,066,055  1,056,334,498       35    149.5      151
#SRR10390685_2.fastq.gz  FASTQ   DNA   7,066,055  1,062,807,718      130    150.4      151

Les fichiers FASTQ contiennent 7 066 055 reads.

Téléchargez le génome de référence de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

# telechargement du genome de ref dans le dir MAPPING
cd /shared/projects/dubii2021/lnegroni/EvaluationM4-M5-main/MAPPING
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.fna.gz

Quelle est la taille de ce génome ?

# regarder a quoi ressemble le fichier => c'est un fichier fasta
less CF_000009045.1_ASM904v1_genomic.fna.gz 
#decompression et utilisation de seqkit
gzip -d GCF_000009045.1_ASM904v1_genomic.fna.gz
module load seqkit
srun seqkit stats --threads 1 GCF_000009045.1_ASM904v1_genomic.fasta
 

La taille de ce génome est de 4,215,606 paires de bases.

Téléchargez l’annotation de la souche ASM904v1 de Bacillus subtilis disponible à cette adresse

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/045/GCF_000009045.1_ASM904v1/GCF_000009045.1_ASM904v1_genomic.gff.gz

Combien de gènes sont connus pour ce génome ?

# la troisiene colonne du fichier .gff3 contient 3 termes : "gene" ou "cds" ou "mRNA". En supposant qu'il ne faut compter que le ternme "gene":

zcat GCF_000009045.1_ASM904v1_genomic.gff.gz |grep -v "#" | cut -f 3 | grep -c 'gene'

4 536 gènes sont recensés dans le fichier d’annotation.

Contrôle qualité

Lancez l’outil fastqc [2] dédié à l’analyse de la qualité des bases issues d’un séquençage haut-débit

#memo: à lancer du dossier EvaluationM4M5-main:
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_1.fastq.gz -o QC/ -t 8
srun --cpus-per-task 8 fastqc FASTQ/SRR10390685_2.fastq.gz -o QC/ -t 8
# resultats enregistrés dans le dossier CLEANING

La qualité des bases vous paraît-elle satisfaisante ? Pourquoi ?

  • [x ] Oui
  • Non

car le score de qualité est supérieur à 30, la longueur des reads supérieure à 100 et les resultats sont comparables entre R1 et R2 comme le montre les fichiers html pour le (read1) [SRR10390685_1_fastqc.html] et le (read2)[SRR10390685_2_fastqc.html]

Lien vers le rapport MulitQC

Est-ce que les reads déposés ont subi une étape de nettoyage avant d’être déposés ? Pourquoi ?

  • Oui
  • Non

car car les reads sont de differentes longueurs

Quelle est la profondeur de séquençage (calculée par rapport à la taille du génome de référence) ?

#taille du genome de reference  = 1 242 608 pb
# profondeur = (nb de read * taille des reads) /taille du genome
#file                    format  type   num_seqs        sum_len  min_len  avg_len  max_len
#SRR10390685_1.fastq.gz  FASTQ   DNA   7,066,055  1,056,334,498       35    149.5      151
#SRR10390685_2.fastq.gz  FASTQ   DNA   7,066,055  1,062,807,718      130    150.4      151
profondeur <- c(2*7066055 * 150 /  4215606)
profondeur

La profondeur de séquençage est de : 503 X.

Nettoyage des reads

Vous voulez maintenant nettoyer un peu vos lectures. Choisissez les paramètres de fastp [3] qui vous semblent adéquats et justifiez-les.

# avec fastp
# filtres = mean quality > 30 longueur >=100, garder seulement les paires
#
module load fastp
fatp --version #fastp 0.20.0
# la ligne de la mort: fastp prends en entrée ("--in") les deux fichiers (car paired end) non nettoyés et crée en sortie ("--out")les fastq nettoyés dans le dir "CLEANING". la fin de ligne correspond aux paramêtres et la sortie est fastp.json est remplacée par fastp.log
srun --cpus-per-task 8 fastp --in1 FASTQ/SRR10390685_1.fastq.gz --in2 FASTQ/SRR10390685_2.fastq.gz --out1 CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz --out2 CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz --html CLEANING/fastp.html --thread 8 --cut_mean_quality 30 --cut_window_size 8 --length_required 100 --cut_tail &> CLEANING/fastp.log

Les paramètres suivants ont été choisis :

Parametre Valeur Explication
cut_mean_quality 30 reduit les erreurs de lecture à 1 sur 1000
cut_window_size 8 fenetre sur laquelle est calculée Q?
length_required 100 les reads doivent être de plus de 100 pb
cut tail supprime les adapateurs?

Ces paramètres ont permis de conserver 13.554096 M (95.909924%) reads pairés, soit une perte de 4.1% % des reads bruts.

Alignement des reads sur le génome de référence

Maintenant, vous allez aligner ces reads nettoyés sur le génome de référence à l’aide de bwa [4] et samtools [5].

#version de bwa: 0.7.17-r1188
module load bwa

# alignement
srun --cpus-per-task=32 bwa mem GCF_000009045.1_ASM904v1_genomic.fasta ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz ../CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 32 > SRR10390685_on_ASM904.1.sam

# charger samtools 1.10
module load samtools

#transformation en fichier binaire (5Go -> 1.4 Go)
srun --cpus-per-task=8 samtools view --threads 8 SRR10390685_on_ASM904.1.sam -b > SRR10390685_on_ASM904.1.bam

#tri et indexation
srun samtools sort SRR10390685_on_ASM904.1.bam -o SRR10390685_on_ASM904.1.sort.bam
srun samtools index SRR10390685_on_ASM904.1.sort.bam

# suppression des fichiers intermediaires via Jupyter Hub
'''les differentes étapes auraient pu être résumées par:
srun --cpus-per-task=34 bwa mem GCF_000009045.1_ASM904v1_genomic.fasta ../CLEANING/SRR10390685_1.cleaned_filtered.fastq.gz ../CLEANING/SRR10390685_2.cleaned_filtered.fastq.gz -t 32 | srun samtools view -b - | srun samtools sort - > SRR10390685_on_ASM904.1.sort.bam
'''

Combien de reads ne sont pas mappés ?

#statistics from samtools
srun samtools idxstats SRR10390685_on_ASM904.1.sort.bam > SRR10390685_on_ASM904.1.sort.bam.idxstats
srun samtools flagstat SRR10390685_on_ASM904.1.sort.bam > SRR10390685_on_ASM904.1.sort.bam.flagstat
# nv de reads non mappés dans le fichier .idxstats
#13571369 + 0 in total (QC-passed reads + QC-failed reads)
#12826829 + 0 mapped (94.51% : N/A)
# soit 744 540 non mappés (5.5%)

744 540 (5.5%) reads ne sont pas mappés.

Croisement de données

  • Calculez le nombre de reads qui chevauchent avec au moins 50% de leur longueur le gène trmNF grâce à l’outil bedtools [6]:
# extraire l'annotation de trmNF
grep   "trmNF" GCF_000009045.1_ASM904v1_genomic.gff | awk '$3=="gene"'> trmNF.gff

# pour voir a quoi ca ressemble:
nano trmNF.gff

#recuperer la sequence genomique correspondant à trmNF
module load bedtools
srun bedtools getfasta -fi GCF_000009045.1_ASM904v1_genomic.fasta  -bed trmNF.gff > trmNF.fasta

#recuperer les reads -f 0.5 pour le chevaucheent à 50% min

srun bedtools intersect -a SRR10390685_on_ASM904.1.sort.bam -b trmNF.gff -sorted -f 0.5 >  trmNF.bam

# tri (sort) index et stat avec samtools
srun samtools sort trmNF.bam -o trmNF.sort.bam
srun samtools index trmNF.sort.bam
srun samtools idxstats trmNF.sort.bam > trmNF.sort.bam.idxstats #pas de nom pour les col de ce fichier...
srun samtools flagstat trmNF.sort.bam > trmNF.sort.bam.flagstat

2801 reads chevauchent le gène d’intérêt.

Visualisation

# pour la visualisation
samtools faidx GCF_000009045.1_ASM904v1_genomic.fasta
samptools faidx SRR10390685_on_ASM904.1.sort.bam
samptools faidx trmNF.sort.bam
'''
download des fichiers 
GCF_000009045.1_ASM904v1_genomic.fasta
GCF_000009045.1_ASM904v1_genomic.fasta.ai
SRR10390685_on_ASM904.1.sort.bam
SRR10390685_on_ASM904.1.sort.bam.ai
trmNF.sort.bam
trmNF.sort.bam.ai

installation de IGV sous win10
loading du fichier genomique (fasta) , son annotation (gff) des bam du jeux complet et de trmNF
decouverte du soft et export image de la zone correpondant à trmNF
'''

Utilisation de IGV local pour la visualisation: igv_snapshot

References

1. toolkit NS. NCBI SRA toolkit. NCBI, GitHub repository. 2019.
2. Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
3. Zhou Y, Chen Y, Chen S, Gu J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34:i884–90. doi:10.1093/bioinformatics/bty560.
4. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:13033997. 2013.
5. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.
6. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
 

A work by Migale Bioinformatics Facility

https://migale.inrae.fr

Our two affiliations to cite us:

Université Paris-Saclay, INRAE, MaIAGE, 78350, Jouy-en-Josas, France

Université Paris-Saclay, INRAE, BioinfOmics, MIGALE bioinformatics facility, 78350, Jouy-en-Josas, France